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Abstract 

We present source code for the computer algebra system Mathematica that analyzes the motion 
of the Pioneer spacecraft using the public available ephemeris data from JPL's website. Within 15 
04 ' minutes, the reader can verify that the Pioneer anomalous acceleration ap (1) exists in the order 

^ ' of magnitude of cHq, (2) is not due to mismodeling of gravitational attraction, solar pressure or 

spacecraft attitude maneuvers. The simple code of about 100 lines may easily be extended by the 
reader to include further tests. Due to the limitations of our approach, we do not know (1) whether 
the unknown raw data were correctly processed to generate the trajectory files (2) how the apparent 
mismatch of ephemerides before 1990 had occurred. 
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1 Introduction 



The Pioneer anomaly ap w —8.7 x 10~^'^ins~^ [1, 2] has become one of the major challenges for theoretical 
physicists, and further efforts to investigate its origin are underway [3, 4, 5, 6]. The approximate coincidence 
of Up with the speed of light divided by the age of the universe has raised the question if this effect is new 
physics or a general failure of Newton's law of gravitation [2, 7, 8, 9]. A couple of months ago, one of us 
heard two astrophysicists talking about the anomaly. One of them, probably tired of hearing about new 
^ ■ trouble besides DM and DE, concluded: 'well, I still don't believe it...'. 

5^ , Though the analysis of the Pioneer data was done by two independent groups and published in peer- 

reviewed journals, it is reality that convincing scientists needs time. Here we do not present any new results 
that help to understand the anomaly, and our approach cannot compete with the detailedness of the 
expert's analysis [2, 7] of the non-public raw data. From a point of view of general scientific methodology, 
we find it however desirable that important results of fundamental physics that require extensive numerical 
treatment can be repeated by a broad public of non-expert scientists. The preliminary analysis and the 
code given below is indeed intended for those who like to get their own opinion in brief. Furthermore, 
minor modifications allow to test some alternative explanations the reader eventually may have in mind. A 
quick description for getting started is found in section 6.1. Though we cannot give a detailed description 
of the program, some clarifying comments are included in the quite self-explaining code (see appendix). 

2 Methods 
2.1 Limitations 

Our analysis is based on the reliability of JPL's ephemeris files. As far as Pioneer 10 and 11 is concerned, 
they contain an explicit warning: 
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'This trajectory is suitable for general historical purposes, but should be used cautiously for 
high precision or tracking data applications. This is due to potential dynamical mismatches 
between the Pioneer era models (DE-118, Lieske's E3 satellite theory of JUP035, SAT050, 
etc.) and the current modern solutions used by Horizons. For example, if the Pioneer 10 (11) 
solutions used here estimated planet or satellite ephemeris corrections at the time. However, the 
transformation from the original DE-118 planetary ephemeris coordinate system to the modern 
frame IS computed by Horizons. (...) The circumstances pertaining to the regeneration of the 
spacecraft trajectory source files are not well known.' 

General remarks on the limitations of the accuracy of JPL ephemerides are found in chap. 19 of the 
HORIZONS manual [10]. 

2.2 General method 

We used the computer algebra system Moihernatica to calculate spacecraft trajectories from known initial 
conditions and from the positions of gravitating planetary bodies. The built-in-function NDSolve with an 
explicit Runge-Kutta method was used to integrate the equations of motion 

^'T) = -G E ^f^^ m = r-,: im = (1) 

while the sum is taken over all relevant bodies. Instead of calculating gravitational forces using the grav- 
itational constant G and masses, the much more accurate Keplerian constants [10], p. 47, for the sun 
and respective planets were implemented. For simplicity. Mercury, Venus, Earth, Mars and the asteroid 
belt masses were assumed to stay at sun's barycenter. Thus, for our approximate analysis ephemerides of 
the outer planets and the sun were sufficient. Since the same 1/r^-law is obeyed, radiation pressure was 
implemented by slightly diminishing the sun's effective mass (see code below for details)"^. We estimated 
an effective surface^ of the spacecraft of 5.9m^ and an albedo of 0.7. 

2.3 Data acquisition 



BODIES OBHITS EPHEMERIDES TOOLS PHYSICAL DATA DISCOVERY FAQ 



HORIZONS Web-Interface 



This tool provides a web-based limited Merisce to JPL's HORIZONS system which can be used to generate ephemerides for 
solar-system bodies. Full access to HORIZONS features is available via the primary telnet interface. HORIZONS system news 
shows recent changes and improvements. A web-interface tutorial is available to assist new users. 



Ciii ieiit Settiiiys 

Ephemeris Type [change] 
Target Body [change] 
Coordinate Origin [change] 
Time Span [change] 
Table Settings [change] 
Display/Output [change] 



VECTORS 

Pioneer 10 Spacecifift [-23] 

Solar System Baiyceiiter (SSBI [500(@0] 

Start=1085-(I1-01, Stop=2003-014»1 , Step= l d 

output units=KM-S; quantities code=2; CSV format=YES; object page= 

download/save (plain text file) 



NO 



Figure 1: Screenshot of HORIZONS site 
^For large distances, the antenna can be assumed to be directed to the sun. 

^Calculated from the diameter of the antenna as given on p. 2 [2]; The value given on p. 28 is different. 
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All the data needed were downloaded from JPL's ephemeris site HORIZONS which contains spacecraft 
trajectories and planetary orbits. With our code, the reader can investigate any period of interest. For 
our analysis, we downloaded the orbits of the sun, of Jupiter, Saturn, Uranus and Neptune^ and the 
trajectories of Pioneer 10 and 11. We chose cartesian coordinates (setting VECTORS) with the solar 
system barycenter as origin. A time step of 1 day was sufficient for our purposes^, though the code works 
with smaller steps without changes. For the numerical treatment, data were 3-D-spline interpolated. Due 
to lack of information, we did not do any maneuver modeling. 

2.4 Modeling 

The beginning of our analysis in 1987 coincides with the analysis carried out by [2], see also [8]. The last 
signal from Pioneer 10 was received in 2002, while the RTG of Pioneer 11 became inoperable in Sept 1995. 
No reasonable analysis can be done after those dates; we do not even know if the very last signals have 
been included for the generation of the trajectory files. The time spans for the downloaded files can be 
chosen from 1985-2003 anyway, since start and end can be chosen separately in simulatespacecraft[\. No 
spacecraft data is available before the respective starts in 1972 and 1973. 

3 Results - a preliminary analysis 

We compare predicted (simulated with conventional gravity) and observed (HORIZONS) values for posi- 
tion, velocity, and acceleration (fig. 2 and 3). All quantities show a significant deviation. 

3.1 Comparison with position 

Here we consider the predicted r{t) of our simulation with the position r^f(t) in the ephemeris file^. 
Velocities and accelerations are obtained by numerical differentiation. Though Pioneer 10 and 11 do not 
move precisely in the direction away from the solar system barycenter, for simplicity the radial components 
only are analyzed. 

An anomalous acceleration is clearly visible for both Pioneer 10 and 11 (bottom graphs). While 
the median (minimizes absolute deviation) is —7.60 x 10~^'^ms~^ and —8.29 x 10~^'^ms~^ respectively, a 
quadratic fit to the acceleration data yields —7.25 x 10~^°ms~^ and —7.05 x 10~^°ms~^ (Pioneer 10 and 
11). 

The anomalous velocity decrease is shown in the middle of fig. (2) and (3). A simple estimate a = ^ 
from the velocity data yields —10.34 x 10~^°ms~^ for Pioneer 10 and —8.72 x 10~^°ms~^ for Pioneer 11. 

The deviation in radial position is shown on top of fig.( 2) and (3). In this case, the estimate 

a = yields -13.7 x IQ-^^ms^^ and -9.52 x W-^"ms-^. 

Of course, different time spans yield slightly different values, as the reader may easily verify. All the 
above values are calculated automatically by our program. 

^Other bodies may easily be included; download the respective ephemeris, add the name to variable planetnames, and 
set the variable plmax to the desired value. 

''Flyby modeling would require much smaller steps and quadrupole moments of planets. 
^While running simulatespacecraft[...], set option vflag to 0. 
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Figure 2: Anomalous deviation (r), velocity (v) and acceleration (a) of Pioneer 10 from 1987-2002. 
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Figure 3: Anomalous deviation (r), velocity (v) and acceleration (a) of Pioneerll from 1987-1995. 
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3.2 Comparison with velocities 

The derivatives v{r) of our predicted r(t) from the simulation are here compared to the velocities vnit) 
from HORIZONS^. Calculating a from Av in this case yields —13.2 x 10~^°ms^^ for Pioneer 10 and 
—4.56 X W~^°ms~^ for Pioneer 11. The anomalous acceleration functions (see bottom graphs of fig. 2 and 3) 
do practically not change. Run simulatespacecraft[-l,...,l] and generateplots[-l] to get the respective plot. 

4 Discussion 

Though the deviation of the observed quantities from the predicted ones are clearly visible and confirm 
the order of magnitude of Up w cHq, there are a couple of results we cannot understand yet. 

First of all, there is a big jump in the data on January 1st, 1990. Before that date, we cannot verify 
the known anomalous acceleration at all, there seems to be an agreement with the predicted trajectories. 
This must be due to a systematic error in the ephemeris programs, i.e. a mismatch between data used 
then and now, as mentioned above in the header of the JPL ephemeris file. 

Further jumps in the acceleration of Pioneer 10, though of much smaller amount, occurred on January 
1st, 1993 and in June 1990. The acceleration remains in the range of Up, however. While these jumps occur 
from a more or less constant level to another, there are a couple of isolated disturbances, see bottom of 
fig. (2) and (3). Most likely those disturbances are due to spacecraft maneuvers we did not model at all. 
Taking the median, those disturbances are practically taken out from the analysis. One should keep in 
mind however that even the estimates from Ar and At; (which were afi'ected by maneuvers) yielded Op 
in the correct order of magnitude. Thus in no case the Pioneer anomaly can be an artifact of maneuver 
mismodeling. 

Looking at the acceleration plots at a high resolution we did not show here, a sinewave disturbance 
with amplitude of about lO^^^m.s^^ appears. The period however varies continuously from about 12 days 
(Pioneer 11. around 1987) to 55 days (Pioneer 10, around 1997) and therefore cannot be attributed to a 
mismodeling of lunar ephemerides or a neglection of tidal effects of the receiver station. It is clearly not a 
manifestation of the annual and diurnal signal published in [2] and must be due to some other systematic 
error. 

5 Conclusions 

Despite the limitations of JPL's ephemeris data, we could verify the anomalous acceleration of Pioneer 
10 and 11, i.e. Up cannot be due to a mismodeling of gravitational attraction, maneuvers or radiation 
pressure. With our code, the reader can easily verify or falsify further hypotheses on the origin of the 
anomaly. To account for the effect of dust, Kuiper or asteroid belt masses, dark energy etc. it suffices 
to run the program with different parameters or to add minor changes to the code. We hope that this 
hands-on demonstration will develop further the interest and critical acceptance of this outstanding effect. 

Acknowledgement. Though we are grateful for any comments, please understand that we cannot 
guarantee functionality or give further support for getting this program to run on your computer. 
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6 Appendix: data preparation and source code 

6.1 Step-by step procedure in 15 minutes 

1. Create your directory 'pioneer' and copy all of the following files in there 

2. Goto NASA's ephemeris site HORIZONS: http://ssd.jpl.nasa.gov/horizons.cgi 

3. Change ephemeris type to VECTORS (vector table) 

4. Set Coordinate origin to [500@0], ecliptic and mean equinox of reference epoch ICRF/ J. 2000.0 

5. Set time span from 1985-01-01 to 2003-01-01, step 1 day 

6. Chose table settings km, km/s, quantities code = 2 (two-state vector x, y, z, vx, vy, vz), CSV format 
= YES, object page = NO. 

7. Display/Output: download/save as plain text file. 

8. Chose target body: sun (sol) 

9. Generate ephemeris and save as sun-85-03.txt. See also screenshot fig. 1. 

10. Proceed in the same manner with Jupiter, Saturn, Uranus, Neptune, Pioneer 10 and 11 as target 
bodies. 

11. Cross-check if the file ending for the Pioneers corresponds to the last two entries of 'spans' (line 2 
of the code). 
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12. Chose 'barycenter' where appropriate (otherwise you'll miss Jupiter's satellites) and save as neptune- 
85-03.txt etc. 

13. Type source from appendix or download it from www.alexander-unzicker.de/pioneer.txt'^ 

14. Change in the first line the path to your pioneer directory 

15. Open a Mathematica *.nb file and run the following commands (see also file end of source): 

SetDirectory["c:\\yourpioneerdirectory"]; (* insert a \\ for any \ in the path*) 

<< pioneer.txt; 

planetini[5,"-85-03"]; 

simulatcspacocraft[-l, 46800. 5(*day start 5.1.87*), 8.75(*insert years*), (*option vflag*)]; 
genererateplots[-l] ; 

16. Proceed likewise with simulatespacecraft[-2, 46798.5(* 3.1.87*), 15.08,0] for Pioneer 10. 
6.2 Source code 

^^^^^^^^^^^^^^^^ gj^^l^g^j^ so"t"t '^^'^^'^^'^^'^^'^^'^^'^^'^^'^^'^^'^^'^^'^^^^'^^^ 

planetnames = {"sun", (*"mcrkur", " vcnus"," earth", "mond"," mars", *)" jupiter", "saturn","uranus", 
" neptune" ," pluto" ," pioneerlO" ," pioneer 1 1" } ; 

spans={"-85-03" ,"-85-03"}; (*fileendings of spacecraft emphemerides*) 
mjd=2400000; (** mean julianian day correction**) 
day=3600*24; (* SI unit is seconds *) 

^ ******** ****** constants for radiation pressure modeling *******) 

CO = 299792458; (*spced of light*) AU = 1.496*10^(11); (* astronomic unit *) 

albedo=0.7; solarK=1367*AU'2/cc*(l-(-albedo); 

surfaceP={5.9,5.9}; (* effective surfaces for radiation pressure in m'2 *) 
massP={258,259}; (* spacecraft masses*) 

(*Instead of masses, the much more accurate Kepler constants are used, see HORIZONS manual p. 47 ****) 
KKs=10^9{132712440017.98698,126712767.857796, 37940626.061137281,5794549.00707, 

6836534.0638792608,981.6008877}; (*kcplcr's constants**) 

(** aU masses inside the asteroid belt (3 10" (21) kg ) added to sun (KKs[[l]]) **) 

KKs[[l]]+=10'9*Apply[Plus, {22032.080486417923, 324858.59882645978,398600.43289693922, 

4902.8005821477636, 42828.314258067119, 200 (*asteroid belt estimate*)}]; 

(** for plots: display year numbers at axes instead of days *) 

subti=4; ti87 =Transpose[{Table[46796.5-|-365.25i/subti,{i,0,16subti}], 

Tablc[If[IntegerQ[i/subti]==True,1987+i/subti,""],{i,0,16subti}]}]; 

ti88 =Transpose[{Table[46796.5+365.25i/subti,{i,0,16 subti}], 

Table[If[IntegerQ[i/subti/2]==True,1987+i/subti,""],{i,0,16subti}]}]; 

SCplots=Table[{0,0,0}, {i,10}]; (* plot variable*) 

SDefaultFont = {"Arial", 8}; 

readplanet[name_,ending_]:=Block[{qwc, wer}, filename=planetnames[[name]] < >ending < >".txt"; 
If[FilcInformation[filename]=={},Print["Missing ephemeris file. Stop."];Break[]]; 
10=3; (* interpolation order of splines *) 
qwe = Import [filename, "CSV"]; 

cutl = Position[qwc, "$$S0E"][[1, 1]]; cut2 = Position[qwe, "$$E0E"][[1, 1]]; 
wer = Take [Drop [qwe, cutl], cut2 - cutl - 1]; 

xyz = 1000*Transpose[Take[Transpose[wer], {3, 5}]]; (*** km - m factor ***) 
'^Do not paste and copy from LaTeX source code, this will create error messages. 
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vxyz = 1000*Transposc[Takc[Transpose[wer], {6, 8}]]; 
time = Flatten[Transpose[wer][[l]]] - mjd; 

xl=Interpolation[Transpose[{tinie,Transpose[xyz][[l]]}],InterpolationOrder — > 10]; 
yl=Interpolation[Transposc[{timc,Transposc[xyz] [[2]]}],IntcrpolationOrdcr — > 10]; 
zl=Interpolation[Transpose[{time,Transpose[xyz][[3]]}],InterpolationOrder — > 10]; 
(** x,y,z, as function, interpolated***) 

vxl=Interpolation[Transposc[{timc,day*Transposc[vxyz][[l]]}],IntcrpolationOrdcr — > 10]; 
vyl=Interpolation[Transpose[{tiine,day*Transpose[vxyz][[2]]}],InterpolationOrder — > 10]; 
vzl=Interpolation[Transpose[{time,day*Transpose[vxyz][[3]]}],InterpolationOrder — > 10]; 

]; 

(* simple differentiation routine for interpolating functions of different range, step: 1 day*) 
numdiff[tab_,st_(*start time in JD*),en_(*end time in JD*)]:=Block[{tablp,tablm,tab2m, tab2p}, 
tablp = RotatcRiglit[tab, 1]; tablm = RotatcRight[tab, -1]; 

dtab = Flatten[{tab[[2]]-tab[[l]],Drop[Drop[(tablm- tablp)/2, {!}], {-l}],tab[[-l]]-tab[[-2]]}/day]; 
ttime = Table[i, {i, st, en, 1}]; 
dtabtime = Transposc[{ttimc, dtab}]; 
dlist=Interpolation[dtabtime] ;dlist] ; 
planetini [plmax_,plfileending_] : =Block [{ } , 

(***** Reading data of all relevant planets 1 sun 2 jup, 3 sat, 4 Ura, 5 nep 

xx=yy=zz={}; 

Print [" Initialize planets ..."]; 

For[k=l,k < =plmax,k++,rcadplanct[k,plfilccnding];AppcndTo[xx,xl];AppcndTo[yy,yl];AppendTo[zz,zl]]]; 
simulatespacecraft[craftflag_,start_,yrs_,vflag_(* set to 1 for comparing observed velocities*)]:=Block[{}, 

^^^^^^^^^^^ g-^Q^^-^j^ j^^^ "V'^'luiOS '^^'^^'^^'^^'^^'^^'^^'^^'^^'^^'^^'^^'^^'^^ 

(** -1 for pio 11, -2, for pio 10, anything else runs automatically*) 
fileending=spans[[craftflag]] ; 
sta= start; 

(* radiation pressure modelled as missing solar mass *) 

KKs [ [1] ] - =solarK *surfaceP [ [craftflag] ] /mass? [ [craftflag] ] ; 

end=Ceiling[start+yrs*365.25]+.5; 

readplanct [craftflag, fllccnding]; 

(**** starting values for simulation ***) 

{xO,yO,zO}={xl[t],yl[t],zl[t]} /. t- >start; 

{vxO,vyO,vzO}={vxl[t],vyl[t],vzl[t]} /. t- >start; (* new*) 

Print ["computing trajectory for ", planetnames[ [craftflag]], " ... "]; 

(***** numeric integration of the equations of motion by Runge-Kutta****) 

plmax=Lcngth [zz] ; 

loe = NDSolve[{ 

x"[t]==Apply[Plus,Table[-day^2KKs[[i]](x[t]-xx[[i]][t])/((x[t]-xx[[i]][t])^2+ 

(y[t]-yy[[i]][t])^2+(z[t]-zz[[i]][t])~2)~(3/2), {i, 1, plmax}]], 

y"[t]==Apply[Plus,Table[-day"2KKs[[i]](y[t]-yy[[i]][t])/((x[t]-xx[[i]][t])~2+ 

(y[t]-yy[[i]][t])'2+(z[t]-zz[[i]][t])^2)^(3/2),{i,l, plmax}]], 

z"[t]==Apply[Plus,Table[-day"2 KKs[[i]](z[t]-zz[[i]][t])/((x[t]-xx[[i]][t])^2+ 

(y[t]-yy[[i]][t])'2+(z[t]-zz[[i]][t])'2)'(3/2), {i, 1, plmax}]], 

x[start]==xO,y[start]==yO,z[start]==zO, 

x' [start] ==vxO,y' [start] ==vyO, z' [start] ==vzO}, 

{x[t], y[t],z[t]}, {t, start, end},Method- >"ExplicitRungeKutta"]; 

r[t.] := Sqrt[loe[[l, 1, 2]]^2 + loe[[l, 2, 2]]'2 + loe[[l, 3, 2]]~2]; (* solution r(t) *) 

v[t_]:=D[r[t],t]/day; 

a[t_]:=D[v[t],t]/day; 
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readplanct[craftflag, filccnding] ; (* reading spacecraft data again*) 

radius[t_]:= Sqrt[xl[t]~2 + yl[t]'2 + zl[t]'2]; (* observed r(t) from Horizons*) 

vradial[t_]:= (vxl[t]*xl[t] + vyl[t]*yl[t] + vzl[t]*zl[t])/day/radius[t]; (*radial component*) 

vobserv=Intcrpolation[Transposc[{Table[i,{i, start, end, l}],Table[vradial[t],{t, start, end, 1}]}]]; 

ranom[t_]:=radius[t]-r[t]; (* anomaly of position **) 

(**** vflag=l compares to observed velocity, vflag=0, to position***) 

If[vflag==l, vvv=numdiff[Tablc[r[t], {t, start, end, 1}], start, end]; 

Clear[vanom];vanom[t_]:=vobserv[t]-vvv[t],vanom=numdiff[Table[ranom[t], {t, start, end, l}],start,end]]; 
aanom=numdiff [Table [vanom[t], {t, start, end, l}],start,end]; 

ap=((vanom[t]/.t — >end)-(vanom[t]/.t — >start))/(end-start)/day; (* ap estimate from velocities *) 
ap2=2((ranom[t]/.t - >end)-(ranom[t]/.t - >start))/(end-start)'2/day'2; (* ap from position ***) 
Print ["anomalous acceleration (from v and r): ", ap, " ", ap2]; 
atab=Tablc[aanom[t], {t, sta, end}]; 
median=Median[atab] ; 

abwei=FindMinimum[Sum[Abs[x - atab[[i]]], {i, Lengtli[atab]}], {x, 0}][[2, 1, 2]]; 
abwci2=FindMinimum[Sqrt[Sum[Abs[x - atab[[i]]] '2, {i, Lengtli[atab]}]], {x, 0}][[2, 1, 2]]; 
Print ["Median, absolute, quadratic deviation minimized: ", median, " ", abwei, " ", abwei2]; 
KKs[[l]]+=solarK*surfaceP[[craftflag]]/massP[[craftflag]]]; (*take out again radiation from solar mass*) 
generateplots[body_]:=Block[{},AO=sta+400;(* clean plot: axes origin shifted days in the future *) 
tic=If [body==- 1 ,ti87,ti88] ; 

SCplots[[body,l]]={Plot[radius[t]-r[t],{t, sta, end}, Ticks— >{tic, Automatic}, 

AxesOrigin- >{A0, 0},AxcsLabel - >{"", "r"},PlotRangc- >A11]]}; 
SCplots[[body,2]]={Plot[vanom[t],{t, sta, end}, Ticks— >{tic, Automatic}, 
AxesOrigin- >{A0, 0},AxesLabel - >{"", "v"},PlotRange- >A11]}; 
SCplots[[body,3]]={PIot[aanom[t],{t, sta, end}, Ticks— >{tic. Automatic}, 
AxesOrigin- >{A0, 0},AxesLabel - >{"", "a"}]}]; 
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